Detailing the Regime Shifts in Maine Coastal Current Behavior
Published
March 14, 2025
Code
{library(sf) library(fvcom) library(tidyverse)library(gmRi)library(patchwork)library(rnaturalearth)library(showtext)library(ncdf4)# Cyclic color palettes in scico# From: https://www.fabiocrameri.ch/colourmaps/library(scico)library(legendry)}# namespace conflictsconflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")# Set the themetheme_set(theme_classic() +map_theme())# Project pathslob_ecol_path <-cs_path("mills", "Projects/Lobster ECOL")fvcom_path <-cs_path("res", "FVCOM/Lobster-ECOL")poly_paths <-cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")# Shapefilesnew_england <-ne_states("united states of america", returnclass ="sf") %>%filter(postal %in%c("VT", "ME", "RI", "MA", "CT", "NH", "NY", "MD", "VA", "NJ", "DE", "NC", "PA", "WV"))canada <-ne_states("canada", returnclass ="sf")# # # Support functions for FVCOMsource(here::here("R/FVCOM_Support.R"))
Code
# Use GMRI styleuse_gmri_style_rmd()
Maine Coastal Current Regime Shift Evaluation
The Maine Coastal Current can be split into a western branch (the western maine coastal current WMCC) & an eastern branch (the eastern Maine coastal current). The two branches are separated by an area off the coast of Penobscot Bay. This is a junction where a portion of the MCC will at times deflect to the South, a seasonal behavior that is influenced by river discharge.
Code
# Project pathfvcom_processed_path <-cs_path("mills", "Projects/Lobster ECOL/FVCOM_processed/MaineCoastalCurrent")# Load the daily surface currents (u,v)daily_mcc_sc <-read_csv(here::here("local_data/MCC_PC_timeseries/MCC_daily_surface_velocities.csv"))# Load the PCA Timeseriesmcc_pc_monthly <-read_csv(here::here("local_data/MCC_PC_timeseries/monthly_surface_velocities_PCA_timeseries.csv"))
Code
# Loading FVCOM# Get some file namesfvcom_surfbot_files <-setNames(list.files(fvcom_path, full.names = T, pattern =".nc"),str_remove(list.files(fvcom_path, full.names = F, pattern =".nc"), ".nc"))# Open onegom3_early <-nc_open(fvcom_surfbot_files[[1]])# Get the mesh itself as a simple feature collectiongom3_mesh <-get_mesh_geometry(gom3_early, what ='lonlat')
Code
#Original Areamcc_turnoff_coords <-tribble(~"lon", ~"lat",-69.34, 43.21, # Bottom Left-69.78, 43.62, # Off Popham-67.45, 44.34, # Off Jonesport-67.4, 43.8, # Bottom right-69.34, 43.21, # Bottom left again)# Make it a polygonmcc_turnoff_poly <-st_polygon(list(cbind(mcc_turnoff_coords$lon, mcc_turnoff_coords$lat))) %>%st_sfc(crs =4326) %>%st_as_sf() %>%mutate(area ="Maine Coastal Current Region") st_geometry(mcc_turnoff_poly) <-"geometry"# Project polygonpoly_projected <-st_transform(mcc_turnoff_poly, st_crs(gom3_mesh)) # Turn off s2sf::sf_use_s2(FALSE)# Trim itmcc_studyarea_mesh <-mesh_trim(mesh = gom3_mesh, domain =st_as_sf(poly_projected) )
Long-term Monthly Current Directions
Before proceeding too far, here is what the average flow direction is across the domain of interest to us:
The two principal components we’re using are derived from daily FVCOM surface-layer current vector information. These first two components explain the following proportions of the variance in the current directions:
Code
# Load the PCA we ran:mcc_pca <-readRDS(here::here("local_data/MCC_PC_timeseries/daily_surface_velocities_PCA.rds"))# Summary - Proportion of variancemcc_pca_summ <-summary(mcc_pca)mcc_pca_summ$importance[1:2, 1:2] %>%round(2) %>%as.data.frame() %>%rownames_to_column("Feature") %>% gt::gt()
Feature
PC1
PC2
Standard deviation
21.46
18.97
Proportion of Variance
0.39
0.31
The loadings for each principal component show these patterns when plotted in space:
We can see from these maps that each mesh element contributed a similar amount to each principal component, which East/West velocities contributing more to PC1 & North/South contributing more to PC2.
Further investigation into how these relate to the southward turnoff flow can be found in MCC_Daily_Workup.Qmd
STARS Regime Shifts - {rstars} Implementation
To investigate potential regime shift dynamics, each monthly principal component timeseries will be evaluated using the STARS methodology.
The following results extend the detailed approach of Stirnimann et al. 2019’s {rstars} repository.
Code
# Stirnimann used these values in their paper:# l = 5, 10, 15, 17.5 years, with monthly data# Huber = 1# Subsampling = (l + 1) / 3# Load the function(s)source(here::here("rstars-master","rSTARS.R"))
Regime Shifts in Maine Coastal Current Principal Components
# Summarise the breakpoint locationsshift_points <- mcc_rstars %>%filter(RSI !=0) %>% dplyr::select(time, PC, var, shift_direction)# Plot the breaks over the monthly dataggplot() +geom_vline(data = shift_points,aes(xintercept = time,color = shift_direction),linewidth =1.5) +geom_line(data = mcc_rstars,aes(time, vals, group = PC),linewidth =0.2, alpha =0.5) +geom_line(data = mcc_rstars,aes(time, vals_rollmean, group = PC),linewidth =0.8, alpha =0.75) +scale_color_gmri() +facet_wrap(~var*PC, ncol =1, labeller =label_wrap_gen()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",title ="STARS Regime Shifts, Monthly Maine Coastal Current")
Export
Code
# # Save itwrite_csv( mcc_rstars, here::here("rstars_results/mcc_monthly_shifts.csv"))
Source Code
---title: "Maine Coastal Current RSTARS"description: | Detailing the Regime Shifts in Maine Coastal Current Behaviordate: "Updated on: `r Sys.Date()`"format: html: code-fold: true code-tools: true df-print: kable self-contained: trueexecute: echo: true warning: false message: false fig.align: "center" comment: ""---```{r}{library(sf) library(fvcom) library(tidyverse)library(gmRi)library(patchwork)library(rnaturalearth)library(showtext)library(ncdf4)# Cyclic color palettes in scico# From: https://www.fabiocrameri.ch/colourmaps/library(scico)library(legendry)}# namespace conflictsconflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")# Set the themetheme_set(theme_classic() +map_theme())# Project pathslob_ecol_path <-cs_path("mills", "Projects/Lobster ECOL")fvcom_path <-cs_path("res", "FVCOM/Lobster-ECOL")poly_paths <-cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")# Shapefilesnew_england <-ne_states("united states of america", returnclass ="sf") %>%filter(postal %in%c("VT", "ME", "RI", "MA", "CT", "NH", "NY", "MD", "VA", "NJ", "DE", "NC", "PA", "WV"))canada <-ne_states("canada", returnclass ="sf")# # # Support functions for FVCOMsource(here::here("R/FVCOM_Support.R"))``````{r}#| label: style-sheet#| results: asis# Use GMRI styleuse_gmri_style_rmd()``````{r}#| label: fonts-config#| echo: false# Path to the directory containing the font file (replace with your actual path)font_dir <-paste0(system.file("stylesheets", package ="gmRi"), "/GMRI_fonts/Avenir/")# Register the fontfont_add(family ="Avenir",file.path(font_dir, "LTe50342.ttf"),bold =file.path(font_dir, "LTe50340.ttf"),italic =file.path(font_dir, "LTe50343.ttf"),bolditalic =file.path(font_dir, "LTe50347.ttf"))# Load the fontshowtext::showtext_auto()```# Maine Coastal Current Regime Shift EvaluationThe Maine Coastal Current can be split into a western branch (the western maine coastal current WMCC) & an eastern branch (the eastern Maine coastal current). The two branches are separated by an area off the coast of Penobscot Bay. This is a junction where a portion of the MCC will at times deflect to the South, a seasonal behavior that is influenced by river discharge.```{r}# Project pathfvcom_processed_path <-cs_path("mills", "Projects/Lobster ECOL/FVCOM_processed/MaineCoastalCurrent")# Load the daily surface currents (u,v)daily_mcc_sc <-read_csv(here::here("local_data/MCC_PC_timeseries/MCC_daily_surface_velocities.csv"))# Load the PCA Timeseriesmcc_pc_monthly <-read_csv(here::here("local_data/MCC_PC_timeseries/monthly_surface_velocities_PCA_timeseries.csv"))``````{r}#| label: load-daily-fvcom# Loading FVCOM# Get some file namesfvcom_surfbot_files <-setNames(list.files(fvcom_path, full.names = T, pattern =".nc"),str_remove(list.files(fvcom_path, full.names = F, pattern =".nc"), ".nc"))# Open onegom3_early <-nc_open(fvcom_surfbot_files[[1]])# Get the mesh itself as a simple feature collectiongom3_mesh <-get_mesh_geometry(gom3_early, what ='lonlat') ``````{r}#| label: assemble the mcc mesh area for mapping#Original Areamcc_turnoff_coords <-tribble(~"lon", ~"lat",-69.34, 43.21, # Bottom Left-69.78, 43.62, # Off Popham-67.45, 44.34, # Off Jonesport-67.4, 43.8, # Bottom right-69.34, 43.21, # Bottom left again)# Make it a polygonmcc_turnoff_poly <-st_polygon(list(cbind(mcc_turnoff_coords$lon, mcc_turnoff_coords$lat))) %>%st_sfc(crs =4326) %>%st_as_sf() %>%mutate(area ="Maine Coastal Current Region") st_geometry(mcc_turnoff_poly) <-"geometry"# Project polygonpoly_projected <-st_transform(mcc_turnoff_poly, st_crs(gom3_mesh)) # Turn off s2sf::sf_use_s2(FALSE)# Trim itmcc_studyarea_mesh <-mesh_trim(mesh = gom3_mesh, domain =st_as_sf(poly_projected) )```#### Long-term Monthly Current DirectionsBefore proceeding too far, here is what the average flow direction is across the domain of interest to us:```{r}#| label: plot current direction#| fig-height: 10# Quick verification that directions make senseaverage_directions <- daily_mcc_sc %>%mutate(period = lubridate::month(time, label =TRUE)) %>%bind_rows(mutate(daily_mcc_sc, period ="Year Round")) %>%mutate(period =factor(period, levels =c(month.abb, "Year Round"))) %>%group_by(period, lonc, latc) %>%summarise(across(.cols =c(u,v), .fns =~mean(.x, na.rm = T))) %>%mutate(# Radian to degree conversionangle =atan2(v, u) *180/ pi, # Convert from radians to degreesangle =if_else(angle<0, angle+360, angle),speed =sqrt(u^2+ v^2), # Compute speed (magnitude)dx = u / speed *0.05, # Scale x-component for visualizationdy = v / speed *0.05) # Scale y-component for visualization# Map themyr_round_map <- average_directions %>%filter(period =="Year Round") %>%ggplot() +geom_sf(data = new_england) +geom_segment(aes(lonc, latc, xend = lonc + dx, yend = latc + dy, color = angle),arrow =arrow(length =unit(0.075, "cm")),linewidth =0.5) +coord_sf(xlim =c(-70, -67.3), ylim =c(43.2, 44.7), expand = F) +scale_y_continuous(breaks =seq(43, 45, .5)) +scale_x_continuous(breaks =seq(-71, -67, 1)) +facet_wrap(~period, ncol =4) +scale_color_scico(palette ='romaO',breaks =c(45, 90, 135, 180, 225, 270, 315, 360),labels =c("NE", "N", "NW", "W", "SW", "S", "SE", "E"),limits =c(0,360)) +guides(color =guide_colring(reverse = T, start = pi/2)) +theme(legend.position ="bottom", legend.title.position ="top") +labs(x ="Longitude", y ="Latitude",color ="Current Flow")yr_round_map``````{r}#| fig-height: 8# Map the months independentlymonthly_flow <- average_directions %>%filter(period !="Year Round") %>%ggplot() +geom_sf(data = new_england) +geom_segment(aes(lonc, latc, xend = lonc + dx, yend = latc + dy, color = angle),arrow =arrow(length =unit(0.03, "cm")),linewidth =0.25) +coord_sf(xlim =c(-70, -67.3), ylim =c(43.2, 44.7), expand = F) +scale_y_continuous(breaks =seq(43, 45, .5)) +scale_x_continuous(breaks =seq(-71, -67, 1)) +facet_wrap(~period, ncol =3) +scale_color_scico(palette ='romaO',breaks =c(45, 90, 135, 180, 225, 270, 315, 360),labels =c("NE", "N", "NW", "W", "SW", "S", "SE", "E"),limits =c(0,360)) +guides(color =guide_colring(reverse = T, start = pi/2)) +theme(legend.position ="none") +labs(x ="Longitude", y ="Latitude",color ="Current Flow")monthly_flow```### Review the MCC Principal ComponentsThe two principal components we're using are derived from daily FVCOM surface-layer current vector information. These first two components explain the following proportions of the variance in the current directions:```{r}# Load the PCA we ran:mcc_pca <-readRDS(here::here("local_data/MCC_PC_timeseries/daily_surface_velocities_PCA.rds"))# Summary - Proportion of variancemcc_pca_summ <-summary(mcc_pca)mcc_pca_summ$importance[1:2, 1:2] %>%round(2) %>%as.data.frame() %>%rownames_to_column("Feature") %>% gt::gt()```The loadings for each principal component show these patterns when plotted in space:```{r}#| fig.height: 6# Pull + Reshape the Rotations / Loadingsmcc_loadings <-data.frame("PC1"= mcc_pca$rotation[,"PC1"],"PC2"= mcc_pca$rotation[,"PC2"]) %>%rownames_to_column("loc") %>%separate(col ="loc", into =c("var", "elem"), sep ="_") %>%pivot_longer(cols =starts_with("PC"), names_to ="PC", values_to ="PC_rotation") %>%mutate(longname =if_else(var =="u", "Eastward Velocity (u)", "Northward Velocity (v)"))# Map The Loadingsloadings_map <- mcc_studyarea_mesh %>%mutate(elem =as.character(elem)) %>%left_join(mcc_loadings) %>%st_as_sf() %>%ggplot() +geom_sf(data = new_england) +geom_sf(aes(fill = PC_rotation)) +facet_grid(PC~longname, labeller =label_wrap_gen(width =10)) +coord_sf(xlim =c(-70.2, -67.1), ylim =c(43.2, 44.7), expand = F) +scale_y_continuous(breaks =seq(43, 45, .5)) +scale_x_continuous(breaks =seq(-71, -67, 1)) +scale_fill_distiller(palette ="RdBu", limits =c(-0.06, 0.06)) +map_theme(legend.title.position ="top",legend.title =element_text(hjust =0.5),legend.position ="right",legend.direction ="horizontal") +labs(fill ="Principal Component Rotation",title ="MCC Principal Component Loadings/Rotations")# Plot daily with the monthly smoothmcc_pc_ts <- mcc_pc_monthly %>%ggplot() +geom_line(aes(date, vals, color =`Principal Component`),linewidth =0.4, alpha =0.4) +geom_line(aes(date, vals_rollmean, color =`Principal Component`),linetype =1) + rcartocolor::scale_color_carto_d() +labs(title ="MCC Surface Current PCA Timeseries",y ="Principal Component Value",x ="Date",color ="12-month average")loadings_map / mcc_pc_ts```We can see from these maps that each mesh element contributed a similar amount to each principal component, which East/West velocities contributing more to PC1 & North/South contributing more to PC2.Further investigation into how these relate to the southward turnoff flow can be found in `MCC_Daily_Workup.Qmd`# STARS Regime Shifts - {rstars} ImplementationTo investigate potential regime shift dynamics, each monthly principal component timeseries will be evaluated using the STARS methodology.The following results extend the detailed approach of Stirnimann et al. 2019's [{rstars}](https://github.com/LStirnimann/rstars) repository.```{r}# Stirnimann used these values in their paper:# l = 5, 10, 15, 17.5 years, with monthly data# Huber = 1# Subsampling = (l + 1) / 3# Load the function(s)source(here::here("rstars-master","rSTARS.R"))```### Regime Shifts in Maine Coastal Current Principal Components```{r}# Get the for themcc_rstars <- mcc_pc_monthly %>%select(`Principal Component`, time = date, vals) %>%split(.$`Principal Component`) %>%map_dfr(function(x){rstars(data.timeseries =as.data.frame( x[,c("time", "vals")]),l.cutoff =12*7, # Seven yearspValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (12*7+1)/3,returnResults = T) },.id ="PC" ) %>%mutate(var ="Maine Coastal Current",PC =str_replace_all(PC, "PC", "Principal Component "),shift_direction =if_else(RSI>0, "Shift Up", "Shift Down")) %>%group_by(PC) %>%arrange(time) %>%mutate(vals_rollmean = zoo::rollmean(vals, k =12, fill =NA, align ="center")) %>%ungroup()```### Summary Plot```{r}#| fig-height: 10# Summarise the breakpoint locationsshift_points <- mcc_rstars %>%filter(RSI !=0) %>% dplyr::select(time, PC, var, shift_direction)# Plot the breaks over the monthly dataggplot() +geom_vline(data = shift_points,aes(xintercept = time,color = shift_direction),linewidth =1.5) +geom_line(data = mcc_rstars,aes(time, vals, group = PC),linewidth =0.2, alpha =0.5) +geom_line(data = mcc_rstars,aes(time, vals_rollmean, group = PC),linewidth =0.8, alpha =0.75) +scale_color_gmri() +facet_wrap(~var*PC, ncol =1, labeller =label_wrap_gen()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",title ="STARS Regime Shifts, Monthly Maine Coastal Current")```### Export```{r}# # Save itwrite_csv( mcc_rstars, here::here("rstars_results/mcc_monthly_shifts.csv"))```